      subroutine parcelv(ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, kwid, &
         ibar, jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, iphead, iphd2, &
         link, transpos, qom, ico, qpar, pxi, peta, pzta, up, vp, wp, pmu, bxv&
         , byv, bzv, wate, uptilde, vptilde, wptilde, mask, pixx, pixy, pixz, &
         piyy, piyz, pizz, qdnv, number, j12x, j12y, j12z, jmag, jxtilde, jytilde&
         , jztilde, jxs, jys, jzs, pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
         pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1,&
		 uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2,avg)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: nvtx
      integer , intent(in) :: nsp
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer , intent(in) :: ibar
      integer , intent(in) :: jbar
      integer , intent(in) :: kbar
      integer  :: itdim
      real(double)  :: avg
      real(double) , intent(in) :: dt
      real(double) , intent(in) :: transpos
      logical , intent(in) :: anydrift
      logical  :: eandm
      integer , intent(in) :: ijkcell(*)
      integer , intent(in) :: ijkvtx(*)
      integer  :: ijkctmp(*)
      integer  :: iphead(*)
      integer , intent(inout) :: iphd2(*)
      integer , intent(inout) :: link(0:*)
      integer , intent(in) :: ico(0:*)
      real(double) , intent(in) :: qom(*)
      real(double) , intent(in) :: qpar(0:*)
      real(double)  :: pxi(0:*)
      real(double)  :: peta(0:*)
      real(double)  :: pzta(0:*)
      real(double) , intent(in) :: up(0:*)
      real(double) , intent(in) :: vp(0:*)
      real(double) , intent(in) :: wp(0:*)
      real(double) , intent(in) :: pmu(0:*)
      real(double) , intent(in) :: bxv(*)
      real(double) , intent(in) :: byv(*)
      real(double) , intent(in) :: bzv(*)
      real(double)  :: wate(itdim,27)
      real(double) , intent(inout) :: uptilde(*)
      real(double) , intent(inout) :: vptilde(*)
      real(double) , intent(inout) :: wptilde(*)
      real(double) , intent(inout) :: mask(*)
      real(double) , intent(inout) :: pixx(*)
      real(double) , intent(inout) :: pixy(*)
      real(double) , intent(inout) :: pixz(*)
      real(double) , intent(inout) :: piyy(*)
      real(double) , intent(inout) :: piyz(*)
      real(double) , intent(inout) :: pizz(*)
      real(double) , intent(inout) :: qdnv(itdim,*)
      real(double) , intent(inout) :: number(itdim,*)
      real(double) , intent(out) :: j12x(*)
      real(double) , intent(out) :: j12y(*)
      real(double) , intent(out) :: j12z(*)
      real(double) , intent(inout) :: jmag(*)
      real(double) , intent(inout) :: jxtilde(*)
      real(double) , intent(inout) :: jytilde(*)
      real(double) , intent(inout) :: jztilde(*)
      real(double) , intent(inout) :: jxs(itdim,*)
      real(double) , intent(inout) :: jys(itdim,*)
      real(double) , intent(inout) :: jzs(itdim,*)
      real(double) , intent(inout) :: pixysp2(*)
      real(double) , intent(inout) :: piyzsp2(*)
      real(double) , intent(inout) :: pixzsp2(*)
      real(double) , intent(inout) :: pixxsp2(*)
      real(double) , intent(inout) :: piyysp2(*)
      real(double) , intent(inout) :: pizzsp2(*)
      real(double) , intent(inout) :: pixysp1(*)
      real(double) , intent(inout) :: piyzsp1(*)
      real(double) , intent(inout) :: pixzsp1(*)
      real(double) , intent(inout) :: pixxsp1(*)
      real(double) , intent(inout) :: piyysp1(*)
      real(double) , intent(inout) :: pizzsp1(*)
	  real(double) , intent(inout) :: uxsp1(*),uysp1(*),uzsp1(*),densp1(*)
	  real(double) , intent(inout) :: uxsp2(*),uysp2(*),uzsp2(*),densp2(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: ijkvstep
      integer :: ibp2, jbp2, kbp2, l, n, newcell, ijk, is, newcell_nxt, np
      real(double) :: nu, bxp, byp, bzp, qomdt, up0, vp0, wp0, udotb, denom
      logical :: nomore
!-----------------------------------------------
!
!
!
!
!
!
!      a routine to interpolate particle data to the vertices of grid
!     LAST CHANGED  2/21/94  8:06AM
!
      ibp2 = ibar + 2
      jbp2 = jbar + 2
      kbp2 = kbar + 2
!
      call vtxindx (iwid, jwid, kwid, ijkvstep)
!
!
!
      wate(ijkvtx(:nvtx),:8) = 0.0
!
!
      iphd2(ijkcell(:ncells)) = 0
      ijkctmp(:ncells) = ijkcell(:ncells)
!
!
      newcell = ncells
!
!     **************************************
!
!     set accumulators to zero
!
!     ********************************************
!
!dir$ ivdep
!      do 25 n=1,nvtx
!
!      ijk=ijkvtx(n)
!
      jmag(:ibp2*jbp2*kbp2) = 0.0
      j12x(:ibp2*jbp2*kbp2) = 0.0
      j12y(:ibp2*jbp2*kbp2) = 0.0
      j12z(:ibp2*jbp2*kbp2) = 0.0
      jxtilde(:ibp2*jbp2*kbp2) = 0.0
      jytilde(:ibp2*jbp2*kbp2) = 0.0
      jztilde(:ibp2*jbp2*kbp2) = 0.0
      pixx(:ibp2*jbp2*kbp2) = 0.0
      pixy(:ibp2*jbp2*kbp2) = 0.0
      pixz(:ibp2*jbp2*kbp2) = 0.0
      piyy(:ibp2*jbp2*kbp2) = 0.0
      piyz(:ibp2*jbp2*kbp2) = 0.0
      pizz(:ibp2*jbp2*kbp2) = 0.0
!
!c      pixysp2(ijk)=0.0
!c      piyzsp2(ijk)=0.0
      pixysp2(:ibp2*jbp2*kbp2) = pixysp2(:ibp2*jbp2*kbp2)*avg
      piyzsp2(:ibp2*jbp2*kbp2) = piyzsp2(:ibp2*jbp2*kbp2)*avg
      pixzsp2(:ibp2*jbp2*kbp2) = pixzsp2(:ibp2*jbp2*kbp2)*avg
      pixxsp2(:ibp2*jbp2*kbp2) = pixxsp2(:ibp2*jbp2*kbp2)*avg
      piyysp2(:ibp2*jbp2*kbp2) = piyysp2(:ibp2*jbp2*kbp2)*avg
      pizzsp2(:ibp2*jbp2*kbp2) = pizzsp2(:ibp2*jbp2*kbp2)*avg
      pixysp1(:ibp2*jbp2*kbp2) = pixysp1(:ibp2*jbp2*kbp2)*avg
      piyzsp1(:ibp2*jbp2*kbp2) = piyzsp1(:ibp2*jbp2*kbp2)*avg
      pixzsp1(:ibp2*jbp2*kbp2) = pixzsp1(:ibp2*jbp2*kbp2)*avg
      pixxsp1(:ibp2*jbp2*kbp2) = pixxsp1(:ibp2*jbp2*kbp2)*avg
      piyysp1(:ibp2*jbp2*kbp2) = piyysp1(:ibp2*jbp2*kbp2)*avg
      pizzsp1(:ibp2*jbp2*kbp2) = pizzsp1(:ibp2*jbp2*kbp2)*avg

!
!
      qdnv(:ibp2*jbp2*kbp2,:nsp) = 0.0
      number(:ibp2*jbp2*kbp2,:nsp) = 0.0
      jxs(:ibp2*jbp2*kbp2,:nsp) = 0.0
      jys(:ibp2*jbp2*kbp2,:nsp) = 0.0
      jzs(:ibp2*jbp2*kbp2,:nsp) = 0.0
!
!
    1 continue
      nomore = .TRUE.
!
!
      newcell_nxt = 0
      do n = 1, newcell
         if (iphead(ijkctmp(n)) == 0) cycle
         newcell_nxt = newcell_nxt + 1
         ijkctmp(newcell_nxt) = ijkctmp(n)
      end do
!
      newcell = newcell_nxt
!     write(66,*)'parcelv,',newcell,nsp
!
      if (newcell /= 0) then
         nomore = .FALSE.
!
         call watev (newcell, ijkctmp, iphead, itdim, pxi, peta, pzta, wate)
!
!dir$ ivdep
!
         do n = 1, newcell
!
            ijk = ijkctmp(n)
!
            bxp = (wate(ijk,1)*bxv(ijk+iwid)+(wate(ijk,2)*bxv(ijk+iwid+jwid)+(&
               wate(ijk,3)*bxv(ijk+jwid)+(wate(ijk,4)*bxv(ijk)+(wate(ijk,5)*bxv&
               (ijk+iwid+kwid)+(wate(ijk,6)*bxv(ijk+iwid+jwid+kwid)+(wate(ijk,7&
               )*bxv(ijk+jwid+kwid)+wate(ijk,8)*bxv(ijk+kwid))))))))*transpos
!
            byp = (wate(ijk,1)*byv(ijk+iwid)+(wate(ijk,2)*byv(ijk+iwid+jwid)+(&
               wate(ijk,3)*byv(ijk+jwid)+(wate(ijk,4)*byv(ijk)+(wate(ijk,5)*byv&
               (ijk+iwid+kwid)+(wate(ijk,6)*byv(ijk+iwid+jwid+kwid)+(wate(ijk,7&
               )*byv(ijk+jwid+kwid)+wate(ijk,8)*byv(ijk+kwid))))))))*transpos
!
            bzp = (wate(ijk,1)*bzv(ijk+iwid)+(wate(ijk,2)*bzv(ijk+iwid+jwid)+(&
               wate(ijk,3)*bzv(ijk+jwid)+(wate(ijk,4)*bzv(ijk)+(wate(ijk,5)*bzv&
               (ijk+iwid+kwid)+(wate(ijk,6)*bzv(ijk+iwid+jwid+kwid)+(wate(ijk,7&
               )*bzv(ijk+jwid+kwid)+wate(ijk,8)*bzv(ijk+kwid))))))))*transpos
!
            is = ico(iphead(ijk))
            qomdt = 0.5*qom(is)*dt
!
            up0 = up(iphead(ijk))
            vp0 = vp(iphead(ijk))
            wp0 = wp(iphead(ijk))
!
            udotb = up0*bxp + vp0*byp + wp0*bzp
!
            denom = 1./(1. + (bxp**2 + byp**2 + bzp**2)*qomdt**2)
!
            uptilde(ijk) = (up0 + qomdt*(vp0*bzp - wp0*byp + qomdt*(udotb*bxp))&
               )*denom
!
            vptilde(ijk) = (vp0 + qomdt*(wp0*bxp - up0*bzp + qomdt*(udotb*byp))&
               )*denom
!
            wptilde(ijk) = (wp0 + qomdt*(up0*byp - vp0*bxp + qomdt*(udotb*bzp))&
               )*denom
!
!
         end do


!dir$ ivdep
         do n = 1, newcell
!
                  ijk = ijkctmp(n)
                  is = ico(iphead(ijk))
!
!     calculate the charge density at the vertices
!
                  qdnv(ijk+ijkvstep(1:8),is) = qdnv(ijk+ijkvstep(1:8),is) + qpar(&
                     iphead(ijk))*wate(ijk,1)

!
                  jxs(ijk+ijkvstep(1:8),is) = jxs(ijk+ijkvstep(1:8),is) + qpar(&
                     iphead(ijk))*wate(ijk,1:8)*uptilde(ijk)

!
                  jys(ijk+ijkvstep(1:8),is) = jys(ijk+ijkvstep(1:8),is) + qpar(&
                     iphead(ijk))*wate(ijk,1:8)*vptilde(ijk)

!
                  jzs(ijk+ijkvstep(1:8),is) = jzs(ijk+ijkvstep(1:8),is) + qpar(&
                     iphead(ijk))*wate(ijk,1:8)*wptilde(ijk)

!
                  j12x(ijk+ijkvstep(1:8)) = j12x(ijk+ijkvstep(1:8)) + qpar(&
                     iphead(ijk))*wate(ijk,1:8)*up(iphead(ijk))
!
                  j12y(ijk+ijkvstep(1:8)) = j12y(ijk+ijkvstep(1:8)) + qpar(&
                     iphead(ijk))*wate(ijk,1:8)*vp(iphead(ijk))
!
                  j12z(ijk+ijkvstep(1:8)) = j12z(ijk+ijkvstep(1:8)) + qpar(&
                     iphead(ijk))*wate(ijk,1:8)*wp(iphead(ijk))
!
                  number(ijk,is) = number(ijk,is) + 1d0
!
             if (mod(is, 2).eq.1) then
!
                  pixysp1(ijk+ijkvstep(1:8)) = pixysp1(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*uptilde(ijk)*&
                     vptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  piyzsp1(ijk+ijkvstep(1:8)) = piyzsp1(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*vptilde(ijk)*&
                     wptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  pixzsp1(ijk+ijkvstep(1:8)) = pixzsp1(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*uptilde(ijk)*&
                     wptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  pixxsp1(ijk+ijkvstep(1:8)) = pixxsp1(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*uptilde(ijk)*&
                     uptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  piyysp1(ijk+ijkvstep(1:8)) = piyysp1(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*vptilde(ijk)*&
                     vptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  pizzsp1(ijk+ijkvstep(1:8)) = pizzsp1(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*wptilde(ijk)*&
                     wptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!
            else
!c
                  pixysp2(ijk+ijkvstep(1:8)) = pixysp2(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*uptilde(ijk)*&
                     vptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  piyzsp2(ijk+ijkvstep(1:8)) = piyzsp2(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*vptilde(ijk)*&
                     wptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  pixzsp2(ijk+ijkvstep(1:8)) = pixzsp2(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*uptilde(ijk)*&
                     wptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  pixxsp2(ijk+ijkvstep(1:8)) = pixxsp2(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*uptilde(ijk)*&
                     uptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  piyysp2(ijk+ijkvstep(1:8)) = piyysp2(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*vptilde(ijk)*&
                     vptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!c
                  pizzsp2(ijk+ijkvstep(1:8)) = pizzsp2(ijk+ijkvstep(1:8)) + qpar(iphead(ijk))*wptilde(ijk)*&
                     wptilde(ijk)*(1d0-avg)*wate(ijk,1:8)
!
         endif
!
    end do
!
!
!
!     ************************************************
!
         if (anydrift) then
!
            do l = 1, 8
!
!dir$ ivdep
               do n = 1, newcell
!
                  ijk = ijkctmp(n)
!
                  jmag(ijk+ijkvstep(l)) = jmag(ijk+ijkvstep(l)) + pmu(iphead(&
                     ijk))*wate(ijk,l)
!
               end do
            end do
!
         endif

!
!     ************************************************
!
!     calculate the pressure at t+dt/2
!
!     pressure is a cell-centered variable
!     calculated by NGP interpolation
!
!     ***********************************************
!dir$ ivdep
         do n = 1, newcell
!
            ijk = ijkctmp(n)
!
            pixx(ijk) = pixx(ijk) + qpar(iphead(ijk))*uptilde(ijk)*uptilde(ijk)
!c
            pixy(ijk) = pixy(ijk) + qpar(iphead(ijk))*uptilde(ijk)*vptilde(ijk)
!c
            pixz(ijk) = pixz(ijk) + qpar(iphead(ijk))*uptilde(ijk)*wptilde(ijk)
!c
            piyy(ijk) = piyy(ijk) + qpar(iphead(ijk))*vptilde(ijk)*vptilde(ijk)
!c
            piyz(ijk) = piyz(ijk) + qpar(iphead(ijk))*vptilde(ijk)*wptilde(ijk)
!c
            pizz(ijk) = pizz(ijk) + qpar(iphead(ijk))*wptilde(ijk)*wptilde(ijk)
!
         end do
!
!
!
!
         do n = 1, newcell
            ijk = ijkctmp(n)
            np = iphead(ijk)
            if (np <= 0) cycle
            iphead(ijk) = link(np)
            link(np) = iphd2(ijk)
            iphd2(ijk) = np
!
         end do
      endif
!
      if (.not.nomore) go to 1
!
!dir$ ivdep
      do n = 1, ncells
!
         ijk = ijkcell(n)
!
         iphead(ijk) = iphd2(ijk)
!
         iphd2(ijk) = 0
!
      end do
!
      do is = 1, nsp
         do n = 1, nvtx
            ijk = ijkvtx(n)
            jxtilde(ijk) = jxtilde(ijk) + jxs(ijk,is)
            jytilde(ijk) = jytilde(ijk) + jys(ijk,is)
            jztilde(ijk) = jztilde(ijk) + jzs(ijk,is)
         end do
         do n = 1, nvtx
            ijk = ijkvtx(n)
            jxs(ijk,is) = jxs(ijk,is)/(qdnv(ijk,is)+1.E-10)
            jys(ijk,is) = jys(ijk,is)/(qdnv(ijk,is)+1.E-10)
            jzs(ijk,is) = jzs(ijk,is)/(qdnv(ijk,is)+1.E-10)
         end do
      end do
!
!
      return
      end subroutine parcelv
